Cotunneling theory of inelastic STM spin spectroscopy 
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We propose cotunneling as the microscopic mechanism that makes possible inelastic electron 
spectroscopy of magnetic atoms in surfaces for a wide range of systems, including single magnetic 
adatoms, molecules and molecular stacks. We describe electronic transport between the scanning 
tip and the conducting surface through the magnetic system (MS) with a generalized Anderson 
model, without making use of effective spin models. Transport and spin dynamics are described 
with an effective cotunneling Hamiltonian in which the correlations in the magnetic system are 
calculated exactly and the coupling to the electrodes is included up to second order in the tip-MS 
and MS-substrate. In the adequate limit our approach is equivalent to the phenomenological Kondo 
exchange model that successfully describe the experiments . We apply our method to study in detail 
inelastic transport in two systems, stacks of Cobalt Phthalocyanines and a single Mn atom on CU2N. 
Our method accounts both, for the large contribution of the inelastic spin exchange events to the 
conductance and the observed conductance asymmetry. 
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INTRODUCTION 



The combination of two powerful techniques, Inelas- 
tic Electron Tunneling spectroscopy (lETS) and Scan- 
ning Tunneling Microscope (STM) makes it possible to 
probe inelastic excitations with subatomic resolution. 
The STM-IETS technique was first applied to the study 
of vibrational excitations of single molecules on surfaces^ 
and has more recently been used to study spin excita- 
tions of a single and a few magnetic atoms and molecules 
deposited on surfaces.^"— In STM-IETS, electrons tun- 
nel between the tip and the conducting substrate going 
through the magnetic system. As the bias voltage V is 
increased, a new conduction channel opens whenever eV 
is larger than the energy of some internal excitation of 
the atom, which results in a stepwise increase of the dif- 
ferential conductance dl/dV and a peak or dip in the 
cPl/dV'^. Tracing the evolution of the elementary ex- 
citations as a function of an applied magnetic field and 
fitting to effective spin Hamiltonians permits to infer the 
single ion magnetic anisotropy tensor as well as exchange 
coupling between adjacent atoms and moleculesi^"— 

The lETS-STM technique has been applied to a variety 
of magnetic systems weakly coupled to a conducting sub- 
strate. The list includes a single transition metal atom 
(Mn, Fe, Co) deposited on a single monolayer of CU2N on 
Coppeririii^, to chains of up to 10 Mn atoms on the same 
substrate^, to Fe-phthalocyanine (Fe-PC) molecules on 
oxidized Cu^, to stacks of Co-PC molecules on Pb^d?^ to 
Mn-PC on PbO^i and, more recently, a single Fe atom on 
InSb, a semiconducting substrateji^ For all these systems 
it is possible to describe the spin exchange assisted tun- 
neling, which accounts for the coupling between trans- 
port electrons and the localized spins of the magnetic 
atoms or molecules, with Kondo-like Hamiltonians J^"— 
Whereas this approach successfully describes the main 
experimental results, including the differential conduc- 
tance, as well as effects related to current driven spin 



dynamics and/or spin polarized tip, there are questions 
that can not be addressed using effective spin models: 

1. Why the spin assisted inelastic conductance is com- 
parable to the elastic contribution, in contrast with 
the phonon-assisted inelastic contribution? 

2. What is the microscopic origin of the spin exchange 
tunneling? 

3. Why the inelastic conductance is not always sym- 
metric with respect to the inversion of the bias po- 
larity? 

In this work we provide a theoretical framework to 
model the existing STM-IETS experiments that ad- 
dresses these questions. Our starting point is a gener- 
alized multi-orbital/multi-site Anderson model, in which 
the electrons in the localized orbitals of the magnetic sys- 
tem (MS) are hybridized to the itinerant states of the tip 
and the surface. The states of the MS are calculated 
by exact diagonalization of a microscopic Hamiltonian 
that can include Coulomb repulsion, crystal field and 
spin-orbit coupling. Transport and spin dynamics are de- 
scribed by means of an effective cotunneling Hamiltonian 
in which the coupling to the tip and surface is included 
up to second order. This approach works provided that 
the charging energy of the MS is much larger than the 
temperature, applied bias potential and the electrode in- 
duced broadening of the MS levels. Thus, the MS must 
be in the Coulomb Blockade situation, where the charge 
is a good quantum number and current flows due to quan- 
tum charge fluctuations, known as cotunnelingj^i^i 

When applied to a single orbital Anderson model, the 
effective cotunneling Hamiltonian that we obtain is iden- 
tical to the Kondo model obtained through the stan- 
dard Schrieffer- Wolff transformation , ^^1^^ Our method 
can applied to systems with more than one localized 
orbital, necessary to address most experimentally rele- 
vant systems The effective cotunneling Hamiltonian 



describes transitions between the different many-body 
states of the MS induced by their couphng to the itin- 
erant electrons. This permits to calculate the scattering 
rates, both for the dissipative dynamics of the spin exci- 
tations of the MS coupled to the leads and those leading 
to the current. 

The rest of the paper is organized as follows. In Sec. 
II we present the derivation of the effective Hamiltonian 
and the procedure used to calculate the current, leav- 
ing some of the technical details for the appendix. In 
Sec. Ill we apply our approach to the case of a single 
site Anderson model, which permits to test our approach 
against well established results. In Sec. IV we implement 
our approach to model transport through stacks of CoPc 
molecules^'-. For that matter, we describe the CoPc 
stacks by means of a Hubbard model. In Sec. V we study 
the case of a single Mn adatom on a CU2N surface,— 1^ us- 
ing a multi-orbital Anderson model where Coulomb in- 
teraction, crystal field and spin-orbit coupling in the MS 
are included in the Hamiltonian and treated exactly, by 
means of numerical diagonalization. In section VI we 
summarize our main results. 



II. THEORY 
A. Effective Hamiltonian 

We describe a magnetic system weakly coupled to two 
electrodes, denoted as tip (T) and surface {S) without 
loss of generality, using the following Hamiltonian: 

n^nT + ns + nus + Vtun. (i) 

Here Ht + "Hs correspond to the Hamiltonian of the two 
electrodes, Hms the magnetic system and Vtun the tun- 
neling Hamiltonian. We shall consider the two electrodes 
as free electron reservoirs, i.e., 'Ht + 'Hs = Xla^a/a/c" 
where /j^ (fa) is the creation (annihilation) operator of 
a quasiparticle with single particle number a = {k, rj, a}, 
with momentum fc, electrode rj — T, S and spin pro- 
jection in the quantization direction a. In general, the 
central region has a complicated many-body Hamilto- 
nian that includes Coulomb repulsion, spin-orbit cou- 
pling, crystal field terms and so on. The many-body 
eigenstates of HmS: k:'^) have a well defined number of 
electrons q. Only 3 charge sectors, q = qo, + and — , are 
relevant. The qo corresponds to the ground state of the 
MS. The sectors + and — correspond to the MS with 
an extra electron {qq + 1) and an extra hole (go ~ 1) re- 
spectively. The Hamiltonian of the isolated MS can be 
written as: 

Hms = ^Eq^n\q,n){q,n\. (2) 
The tunneling Hamiltonian is given by 

Vtun = J2 ^^rJa'^i + ^-C- = V" + V+, (3) 
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FIG. 1: (Color onlme) Scheme of the cotunneling transport 
process through an almost degenerate multiplet of orbital lev- 
els far from resonant. (a)Elastic transport process without 
change in the magnetic system state, (b) Inelastic excita- 
tion/relaxation process leading to the creation of an electron- 
hole pair in one of the electrodes and no charge transport, 
(c) Inelastic transport process with change of the magnetic 
system state. In cases (b) and (c), the energy of the new 
configuration of the magnetic system changes from Ed ^ E'^. 



where the tunneling of electrons in and out the MS are 
described by V"*" and V~ respectively. Here d; (c?i) are 
the creation (annihilation) operator of an electron in a 
single particle state i = {i, cr} with orbital quantum num- 
bers i and spin a. We assume that single-particle tun- 
neling events are spin conserving and spin independent. 

We will start with the uncoupled Hamiltonian Hq = 
Ht + Hs + Hms- Since Ho commutes with the charge 
operator of the MS, the eigenstates of Ho can be la- 
beled according to the charge q in the central atom. We 
assume that the eigenvalues in the go sector are sepa- 
rated by a large gap from the states in the q — ± sec- 
tors, see Fig. [TJ In particular, the chemical potentials 
of the MS, defined as Hh = Ec^Ne) - EdNe - I) and 
fie = EaiNe + 1) - EaiNe), with EaiNe) the ground 
state energy corresponding to Ng electrons, must satisfy 
/ih \fJ-e ~ fJ-ij\ ^ khT, This corresponds to the 

conditions of deep cotunneling in which the sequential- 
first order transitions are exponentially suppressedi^"— 
In this limit we can use degenerate perturbation theory 
to determine the dynamics of the states in the go sector, 
which we denote with |iV). These states are tensor prod- 
ucts of the electrode ground states and the many body 
states |g, n) of the magnetic system. The tunneling oper- 
ator ^ connects them to states \M±) that are products 
of electrode states with 1 quasiparticle and MS states 
|go ± l,m). Unless otherwise stated, in the rest of the 
paper we label the MS islands states with the shorthand 
notation \n) = \qo,n) and \m±) = |go ± l,m). 

Using degenerate perturbation^L theory we can obtain 
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an effective Hamiltonian for the qo sector where the tun- 
nehng events are included to the lowest order: 
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(4) 



In the calculation of the effective Hamiltonian we are ne- 
glecting the energy variations of the unperturbed states 
inside the go manifold, all taken to be Eq, compared 
to the charging energy. When expanding this operator 
in the basis of the electrode quasiparticles and the MS 
many-body states, we can write the effective Hamiltonian 
for the qo sector as (see Appendix 1X1 for details): 
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are operators that act exclusively on the subspace |goj 
of the neutral MS. Their matrix elements read: 
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where 



7"^'(M',cro-') = {n\d^„\nl+){m+\dl^,\n') (9) 
7r„'(M',o-cr') = {n\d\^,\m^){m-\d,,„\n'). (10) 



Eqs. ([Sl fTOl) constitute the cornerstone of the formal- 
ism. The Hamiltonian 'Hcotun in Eq. ([5|) describes the 
scattering of a quasiparticle from the single particle state 
a' to a in the electrodes together with a transition be- 
tween two many-body states of the MS within the go 
manifold. Three types of elementary processes are de- 
scribed by the effective cotunneling Hamiltonian: elas- 
tic processes in which transport electrons are transferred 
between both electrodes without changes in the central 
region, creation of electron-hole pair in a given electrode 
with the corresponding transition in the central island, 
and inelastic tunneling events. In all of them, it is appar- 
ent from Eqs. ^ and (fTO| that the excitations within 
the go manifold in the MS occurs via virtual transitions 
to the charged manifolds q = — and q = -\-. An scheme 
of each of these processes can be seen in Fig. [1] 

Very much like in the case of effective Kondo models, 
the quasiparticle scattering events can be classified in 
four groups depending on whether they include, or not, 
spin flip and/or electrode transition. In turn, the spin 



conserving events are split in two more groups, depend- 
ing on weather or not they have spin dependent ampli- 
tudes. Because of the spin rotational invariance imposed 
in the tunneling Hamiltonian (U), quasiparticle spin flip 
events imply spin transfer to the MS. Finally, the last 
term in Eq. ([5]) describes a renormalization of the many- 
body levels of the MS and can be re-adsorbed into a new 

Hamiltonian for the central part, Hf^g = 'Hms+X^q Tao) , 
so it will be omitted in the following analysis. 

For a fixed set of initial and final quasiparticle states, 
a, a', the matrices (HKT]) have, at most, the dimension of 
the go manifold. For instance, as discussed in detail in 
Sec. mil in so called Anderson model, when the states 
with g = go in the island are those of an unpaired elec- 
tron, the dimension of the matrices ([7][8]) is 2, correspond- 
ing to the two spin projections of a spin 1/2. As a result, 
the Hamiltonian [S] describes a Kondo coupling between 
the electrode and the spin 1/2 of the MS. 



B. Master equation, transition rates and current 

The procedure described above yields an effective 
Hamiltonian of the MS coupled to the electrodes for 
which the states of the g = ± sectors have been inte- 
grated out. The effective total Hamiltonian of the elec- 
trodes coupled to the go manifold reads: 
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(11) 



This Hamiltonian serves as starting point to calculate 
both current and dynamics of the many body states of 
the MS within the go manifold. The dissipative dynamics 
of the n states in the MS is induced by the coupling to the 
electrodes as described by Hcotun- The master equation 
for the populations of the MS states, P„, is given by 



dPn 

dt 



(12) 



where the transition rates Wnn' for the MS to go from 
state n to n! due to quasiparticle scattering in the elec- 
trodes are calculated by applying the Fermi Golden Rule 
with the perturbation given by the tunneling Hamilto- 
nian ([5]). The steady state solutions of this master equa- 
tion depend, in general, on the Hamiltonian parameters, 
the temperature and the bias voltage. At zero bias, the 
steady state solutions are those of thermal equilibrium. 
At finite bias, P„(l^) can depart significantly from equi- 
librium depending on the relative efficiency of the trans- 
port assisted excitations and relaxations)^ 

The rates Wn,n' are the sum of scattering processes in 
which the initial and final electrode and spin quantum 
numbers of the quasiparticle are well defined. 
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(13) 
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An explicit expression for the spin and electrode depen- 
dent scattering rate " is given in the Appendix El 
Eq. (jBip . The expression involves a convolution over the 
energy dependent density of states and effective cotun- 
neling rates. A simpler expression is obtained by doing a 
number of approximations^SiSS, as explained in Appendix 
IB|) . First, we assume that the electrodes have a flat den- 
sity of states within a bandwidth larger than all relevant 
energy scales in the problem: temperature, bias and the 
excitations energies of the MS within the 50 manifold. 
Second, we neglect the energy dependence of the hop- 
ping matrix elements V-kr\,i = yr\,i^ These approxima- 
tions are justified in lETS experiments where the tem- 
perature is at most a few Kelvins and the applied bias is 
bellow 50mV. If we introduce the excitation energy as- 
sociated to the transition between n' and n states in the 
go manifold, A„„' — En — En' and we define the average 
energy e^^, — 1/2 (/x^ + /i^/ -I- A„„/), the transition rates 
WllZ obtained in Appendix [B] can be expressed as 
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where Q[yi) 



and p^CT are the spin and elec- 
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trode resolved density of states. The many-body matrix 
elements S^^',** " (t) are given by 



kria,k' 7]' cr' k' r}' a' . krjcr 



,(15) 



where k = k{e), i.e., the quasiparticle energy that appear 
in the denominators are replaced by the corresponding 
bias-dependent average energy e^^, . 

In this context, the current is given byi^i^ 



71, n' 



(16) 



where e is the (negative) electron charge. This equa- 
tion has a physically transparent meaning: the current 
is proportional to the transition rates of quasiparticles 
changing electrode. These rates involve transitions of the 
MS from the state n, which is occupied with probability 
PniV), to state n', including elastic events n = n' . 

Our convention for the applied bias is such that eV — 
IJ'S ~ Mr (electrons move from tip to surface for a positive 
applied bias). The bias implies a small charge accumula- 
tion both in the tip and the surface which in turn involves 
a shift of their chemical potentials with respect to their 
equilibrium value, denoted by Ep- Without loss of gener- 
ality we can write 113 — Ep +xeV, iip = Ep + (x— l)eV, 
where x is an undetermined parameter that relates the 
bias voltage to the shift of the chemical potential in each 
electrode. Given the fact that the capacitance of the sur- 
face is much larger than that of the tip, it is reasonable 
to take a; = 0. As we show below, this assumption makes 
it possible to account for the conductance asymmetry re- 
ported experimentallyi^ii^ 



In the following, we shall express the differential con- 
ductance in units of 



Go 

50 = —PsPT {Jts + Wts) , 



(17) 



where Gq = 2e^ /h is the quantum of conductance, 
Pr, = pria and Jts and Wts are just the general- 
izations of the (momentum independent) exchange and 
direct coupling respectively that appears in the Anderson 
model, as it will be shown bellow: 

Jts = 2yj''V^*') [(/i^ - Ep)-^ + {Ep - ^i^r'] , (18) 
and 

Wts = l/i''V^'"V2 [(Me - Ep)-' - [Ep - ^,)'i] (19) 

where vj}^^^ is the maximum value of the couplings be- 
tween electrode 77 and the orbitals of the MS. 



C. Summary of the method 

The approach described above can be implemented in 
a wide range of situations following a sequence of well 
defined steps: 

1. Diagonalization of the MS Hamiltonian in the 3 
relevant charge sectors, q = qo — 1, go, 9o + 1, 
providing \q,n) and Eq^n- 

2. Computation of the matrix elements ([7]) and ^ of 
the effective tunneling Hamiltonian operator, which 
requires the calculation of the many body matrix 
elements 7 (Eqs. ([9]) and pO|) ) and the O- matrix 
prefactors. 

3. Calculation of the scattering rates (dH), which de- 
pend on bias, temperature, MS-electrode coupling, 
electrode density of states and MS wave functions. 

4. Finding the non-equilibrium steady state solutions 
Pn{V) of the master equation 

5. Evaluation of the current using Eq. (|16p . 



D. Comparison with other cotunneling theories 

The calculation of cotunneling current has been widely 
studied before, using different methodologies, mainly in 
the context of quantum dot a^^'^^i'^^ and, more recently, 
moleculesi^"— For instance, in Ref . [SJ-jS^ they compute 
the cotunneling scattering rates by truncating the T- 
matrix down to second order in the electrode coupling. 
On the other hand, a more formal and accurate treat- 
ment, valid also in the strong-coupling regime, was in- 
troduced in Ref. [ssl . where the non-equilibrium Keldysh 
Green function formalism was used to study the inelastic 
spectroscopy of single adsorbed molecules. 
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Whereas there current obtained using these different 
methods is the same, our approach permits to derive 
an effective Hamiltonian which, in the adequate hmit, 
is the same than the effective Kondo Hamihonian used 
extensively in previous works J^^— An interesting work 
addressing the relation between multiple-impurity An- 
derson model at half-filling and a Kondo model was pre- 
sented in Ref. 39, where authors proved that a Hubbard 
chain of N impurities coupled in parallel can be described 
with a S = N/2 SU{2) spin Kondo model. In Ref. M, 
authors used the same generalized Schrieffer- Wolff trans- 
formation to relate a singlet-triplet Anderson impurity 
with a spin model close to its quantum phase transition. 

Our approach, based on an effective cotunneling 
Hamiltonian directly obtained from the exact description 
of the magnetic system, provides a microscopic justifica- 
tion of earlier phenomenological works, at the time that 
it keeps the simplicity that allows to calculate the current 
as described above. 



III. SINGLE ORBITAL ANDERSON MODEL 

In this section we revisit the very well known Anderson 
modelSS. for which the MS is a single site Hubbard model: 



(20) 



where Ed is the on-site energy level, U the on-site 
Coulomb repulsion and Uia- = d\^dia. We now derive 
an effective cotunneling Hamiltonian which, as we show 
below, turns out to be identical to the spin 1/2 Kondo 
model by means of a Schrieffer- Wolff transformation . ^^'^^ 
By so doing, we test the validity of our approach and shed 
some light on the origin of the large contribution of the 
inelastic spin assisted tunneling to the conductance. 

The single-site Hubbard Hamiltonian has only 3 pos- 
sible charge states, empty, singly and doubly occupied. 
The singly occupied manifold has two states, | t) and | \) 
with energy Ed- The empty and doubly occupied mani- 
folds have only 1 state each, | ti) with energy 2Ed + U 
for the -|- manifold and |0), with energy for the — man- 
ifold respectively, li Ed + U :^ Ep Ed '> hT the 
ground state has qo — I and classical charge fiuctuations 
are frozen. Hence, the virtual transition operators act- 
ing on the qo = 1 space have dimension two and can be 
expressed as Pauli matrices, acting on the spin space. 

After a straightforward calculation we find the effective 
cotunneling Hamiltonian with 3 contributions. First, the 
famous exchange assisted Kondo tcrm i^^i^^ 
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The second term 7^2 = J2kk' ,7j7j' .era' '^2{kk' ,r]r]' , a) 
in the Hamiltonian corresponds to a direct (spin- 
independent) interaction, also obtained in the Schrieffcr- 
Wolff transformation^^ 
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Notice how in this model, the exchange assisted Jkk' ,riri' 
and the direct tunneling term Wkk',riri' have a common 
origin, namely, virtual charging of the magnetic site. Im- 
portantly, we see how we can have the spin-flip term much 
larger than the direct term. In particular, in the so called 
symmetric case, for which Ed + U — Ep = Ep — Ed, 
the direct term vanishes altogether, due to a cancellation 
between the electron addition and hole addition chan- 
nels. In that situation only the spin- flip assisted tun- 
neling would be possible. Thus, the cotunneling picture 
provides a natural scenario for the large contribution of 
the inelastic contribution to the conductance. Finally, 
a third term H3 =^ ^ ^ T-i^ik, rj, a) is obtained, which 
can be considered as a renormalization of the on-site en- 
ergy level. 



IV. STACKS OF CoPc MOLECULES 

In this section we model the lETS experiments of 
stacks of Cobalt phthalocyanine molecules (CoPc) de- 
posited on Pb(lll)<^iiS CoPc molecules are planar 
molecules with symmetry and a single Cobalt (Co) 
atom at its center, surrounded by four Nitrogen neighbors 
and enclosed by aromatic macrocycles. A single CoPc 
has a ground state with spin S — 1/2, corresponding to 
an unpaired electron presumably in the d^2 orbital of 
Co. In a stack with N + 1 CoPc molecules, the CoPc 
in contact with the Pb surface acts as a dead layer that 
isolates the remaining molecules. 

The stacking seems to be such that Cobalt atoms are 
underneath Nitrogen atoms of the adjacent molecule. 
The lETS results^ii2 of stacks with N active CoPc (N+1 
molecules in total) can be interpreted as if the molecules 
are coupled via an antiferromagnetic coupling, which pre- 
sumably comes from super-exchange between two Cobalt 
coupled to a common Nitrogen. The observed spin-flip 
excitations were successfully described using a Heisen- 
berg model with an antiferromagnetic (AF) coupling 
J ~ 18meV . 

Whereas the Heisenberg model accounts for the ob- 
served excitation energies, it can not account for either 
the transport mechanism or the fact that the conduc- 
tance in this system is very asymmetric. In particular, 
some inelastic steps seen at a given bias polarity are not 



seen when bias sign is reversed. Additional experiments 
where the charge state of the molecular stack was con- 
trolled using the STM tip as a local gate make it nec- 
essary to go beyond spin-only models.^- The observed 
excitation energies could be accounted for using a Hub- 
bard model, rather than a Heisenberg model: 

■i(7 i 

+tY,[Aad^+l<y+\^■c) (25) 
i.cr 

Here Ed stands for the energy of the dz2_j.2 orbital with 
respect to the Fermi energy, that we take at 0, U stands 
for the on-site Coulomb repulsion and t for the Co-Co 
hopping, which actually occurs through the common Ni- 
trogen neighbor. In the strongly insulating limit, U ^ t 
and at half- filling (1 unpaired electron per Cobalt atom), 
the Hubbard model has the same low energy excitation 
spectra than the Heisenberg model with J = Away 
from half-filling, when the molecular stack is charged, the 
mapping to the Heisenberg model is no longer possible 
but still the excitation energies observed experimentally 
are accounted for by the Hubbard model.— Here we fo- 
cus on the half-filling case and we apply our formalism 
to short Hubbard chains with iV = 2, 3, 4 sites. We take 
U = 1.5eV which imposes t = 82meV (J w 18meV), in 
accordance with the experimentally observed valuci^iiS 

A. The dimer 

The eigenvalues and eigenvectors of the Hubbard dimer 
can be found analytically, both for the half-filling sector 
and the two sectors with 1 and 3 electrons. At half fill- 
ing (q = qo) the ground state corresponds to a spin sin- 
glet, S' = 0, while the first excited state corresponds to 
a spin triplet, 5 = 1 with excitation energy J neglecting 
terms of order t'^/U^, see Fig. [2jb) (in agreement with 
the experimental results^ii^) . Referred to half- filling, the 
electron addition and hole addition energies are Ed + U 
and —Ed, respectively. Thus, for the MS to be at half 
filling we must have Ed < and \Ed\ < U. The states of 
the q — i sectors correspond to those of a single electron 
and a single hole, respectively. 

In order to assess the effect of the relative weight of the 
cotunneling mediated by virtual hole and virtual electron 
addition {q± channels), we have calculated dl/dV spectra 
obtained for the Hubbard dimer for three different val- 
ues of the on-site energy Ed: Ed = —U/2, the so called 
electron hole symmetry point, Ed = — 0.35eV for which 
virtual transitions to the q- manifold are favored and 
Ed = — 1.3eV, which favors virtual transitions to the q^ 
manifold, see Fig. [TJc). Whereas the excitation step at 
±eV = J is present in all of them, both the magnitude 
and the bias dependence of the elastic contribution de- 
pends a lot on Ed- At the electron- hole symmetry point 
(EHSP), the elastic conductance is zero, as in the Ander- 
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FIG. 2: (Color online) (a) Scheme of a two sites Hubbard 
dimer connected to two electrodes, (b) Lowest energy levels 
of Hubbard dimer at half filling in terms of the exchange con- 
stant J. dl /dV as a function of applied bias. Here psVs = 5 
and ptVt = 1 and T = 0.4K. 



son model, and the non-monotonic lineshape right above 
the inelastic step is due to the depletion of the occupa- 
tion of the ground state in favor of the excited state, a 
non-equilibrium effect discussed in our previous worki^^ 
Both the elastic and the inelastic contributions increase 
when Ed is taken away from the EHSP. 

When Ed = — 1.3eV, the virtual transition to the 
manifold is dominant and cotunneling is mediated by the 
addition of an electron. As we mentioned in Sec. Ill Bl our 
bias convention is such that positive bias V results in an 
increment of the tip chemical potential with respect to 
the molecules and the surface. Thus, for V > it be- 
comes easier to add an electron to the system, increasing 
the global conductance. For V < 0, instead, the chemi- 
cal potential of the tip is decreased, making it relatively 
harder to charge the dimer with an electron and reduc- 
ing the cotunneling conductance thereby. In the case of 
Ed = —0.35eV the situation is reversed. The virtual tran- 
sition to the q- manifold is dominant, ie., cotunneling is 
mediated by the addition of a hole (or the removal of an 
electron). In this case a positive bias makes it harder 
for the electron to tunnel out of the system, decreasing 
the conductance. Thus, in our calculation the asymme- 
try of the conductance comes from the assumption that 
the bias shifts mostly the tip chemical potential, and not 
the surface, and the fact that one of the two cotunneling 
channels (virtual addition of either an electron or a hole) 
is dominant. Comparing with the experimental results,— 
we infer that the double CoPc molecule system is close to 
the electron addition point. This has been further con- 
firmed by additional experiments by the same group^i^ 
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FIG. 3: (Color online) dl /dV as a function of applied bias for 
N = 3 with on-site energy Ed — — 0.35eV (red-dashed line), 
Ed = -0.75eV (thin-black line) and Ed = -1.2eV (thick-blue 
line), (a) Serial Hubbard trimer with Vs,i = 25i,3 and Vr,! ~ 
5(5i,i. (b) Multiple electrode connected Hubbard trimer with 
with psVs ~ (2, 1, 1) and psVs = (2,3,5). Inset: scheme of 
the lowest energy levels. The other parameters are kept as in 

Fig.m 



FIG. 4: (Color online) dl /dV as a function of applied bias for 
= 4 with on-site energy Ed = — 0.35eV (red-dashed line), 
Ed = -0.75eV (thin-black line) and Ed = -1.2eV (thick- 
blue line), (a) Serial Hubbard tetramer with Vs,i = 35i,4 
and Vr.i = 105i,i. (b) Multiple electrode connected Hub- 
bard tetramer with with psVs = (0, 0, 0.2, 10) and pt Vs = 
(3,0.3,0,0). Inset: scheme of the lowest energy levels. The 
other parameters are kept as in Fig. [T] 



B. The trimer and the tetramer 

We now consider the Hubbard chains with either N = 
3 and A'' = 4 sites and try to model the CoPc molecular 
stacks with 3 and 4 active molecules respectively^ We 
assume that t and U take the same values than before 
and that there is one electron per site in the ground state. 
We label the sites from n = 1 to ti = A^, starting from 
the molecule closest to the tip. For = 3, the ground 
state and first excited state have S — 1/2 and the second 
excited state has S — 3/2 (see Fig. Thus, we expect 
two inelastic transitions, at energies J and 3J/2. For 
the A'' = 4 chain, the ground state has 5 = and the 
two lowest energy excited states, both with 5 = 1, have 
excitation energies 0.7 J and 1.4 J, see Fig.|31 Again, two 
inelastic steps are expected at those energies. 

In Figs, inja) and IDJa) we show the conductance for 
iV = 3 and = 4 respectively assuming that the elec- 
trons can tunnel from the tip to the n = 1 site only 
and from the n = N site to the surface only. As in the 
case of the dimer, we take 3 different values for Ed'- hole 
mediated, electron mediated and EHSP. On top of the 
symmetry trends already discussed for the dimer, we see 
how in the EHSP only the lowest energy transition is 
seen both in the A^ = 3 and N = A cases. This suggests 
that not only the elastic contribution vanishes, as in the 
case of the Anderson model and the Hubbard dimer, but 
also some of the inelastic transitions can be suppressed 
possibly due to the destructive interference between the 
hole and electron channels. 

In the case of Figs. |31[a) andlH^a), where only the sites 
at the end of the chain are coupled to either the tip or 



the surface, the steps are visible for both signs of V , at 
odds with the experimental observationsi^ In an attempt 
to explore a scenario in which the height of the steps 
are only visible at a given polarity, we have considered 
a situation where electrons can tunnel from the tip to 
sites other than n = 1 and from the surface to sites other 
than n = N . By so doing, we can obtain dl/dV curves 
where the steps are depleted for V < (Figs, ^h) and 
mb)). However, we think that a more plausible explana- 
tion would come from a microscopic calculation including 
more than 1 orbital per molecule. It must also be men- 
tioned the broadening of the excitations observed exper- 
imentally is larger than 5.4 fc^T, which indicates than 
neglecting the intrinsic broadening due to the coupling 
to the continuum of states of the electrodes is not fully 
justified4S 



V. MAGNETIC ADATOMS 

We now consider spin lETS through a single Mn 
atom deposited on a CU2N surface. This sys- 
tem has been widely studied experimentally»3-i4!S and 
theoreticallyfi^"— i^i^ in most instances modeling the 
Mn spin with an effective spin model. Here we go beyond 
the spin model picture and we use a multiorbital Ander- 
son Hamiltonian for the 5 d electrons of the Mn+^ ion 
which includes Coulomb interaction, spin-orbit coupling 
and crystal field. Transport occurs via virtual transitions 
to the many-body states with either 4 or 6 d electrons. 
Our approach requires the exact diagonalization of the 
fermionic model in the 3 relevant charge states, with 4, 
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5 and 6 electrons. Below we describe the multi-orbital 
Anderson model, the transport calculation and compare 
with the experimental results in Ref. 0- 

A. Magnetic system Hamiltonian 

Here we describe our model Hamiltonian for the Mn 
ion in the CU2N surface. The purpose of our model is to 
provide a minimal fermionic Hamiltonian that accounts 
for the data, rather than to provide a realistic descrip- 
tion of the Mn ion on the surface. Density functional 
calculations^i^i^ suggest that the Mn adatom transfers 
charge to the CuN surface and creates bonds with its 
neighboring N atoms. As a result, the Mn adatom be- 
comes a Mn^"*" ion that has lost its two 4s electrons. We 
model this system considering only the electrons of 
the Mn, including the electrostatic potential of the neigh- 
boring atoms. The Hamiltonian of the MS can be written 
as: 

'He = Hee + HcF + Hso + Hzeem, (26) 

where Hee is the Coulomb repulsion between the Sd elec- 
trons, HcF is the crystal field Hamiltonian, Hso is the 
spin-orbit Hamiltonian and Hzeem is the Zeeman Hamil- 
tonian associated to an applied magnetic field B. The 
Coulomb matrix elements of the atomic orbitals can be 
expressed in terms of radial integrals which depend on 
the specific form of the approximate wave function and 
an angular part that can be obtained analytically^^ We 
have taken them from a calculation for an isolated ion 
using Gaussian package^ which yields to the unscreened 
on-site Coulomb repulsion U ~ 24eV, in accordance with 
unscreened Hartree-Fock calculations^^. Since screening 
in the real system makes U much smaller than the single 
ion calculation we have downscaled the Coulomb matrix 
elements with an overall dielectric constant of e = 4.7 in 
order to obtain U in the range of 5eVM 

The energy Ed of the d levels before crystal splitting 
is included, is kept as a free parameter in our theory. 
The crystal field term Hcf, is built using a point charge 
model for the first N and Cu neighbors,'*^ whereas an 
effective dielectric constant e' was introduced to account 
for the screening of the bare crystal field and fit the many- 
body spectrum to that of the single ion Hamiltonian^^S 
Fig. El^a) shows the splitting of the five energy levels 
due to the crystal field, together with its dominant orbital 
contribution. Finally, the spin-orbit Hamiltonian reads: 

Hso^^ J2 (m(7|L-S^|mV)dLdw,.' (27) 

where A = 43meV corresponds to the value of the bare 
Mn^+ ion^ and e" is another free parameter in our 
model. 

The Hamiltonian ([26]) corresponding to the 3d^ elec- 
trons was then diagonalized in the space of the 252 pos- 
sible configurations, using the configuration interaction 



(CI) method. Analogously, the eigenvalues and eigenvec- 
tors of the 3(i^ and 3d^ configurations were calculated in 
order to get the transition rates ()B1[) . The condition of 
stable configuration with Ne = 5 electrons require that 
Eg{5) < Eg{4:), Eg{G)- In our case, this bound trans- 
lates into the inequality -24.1cV < Ed < -18.9cV. In 
particular, we choose Ed in the middle of this energy win- 
dow and, as it will be shown in next section, results do 
not change significantly with Ed- 



B. Mn energy spectra 

According to first Hund's rule, we expect that the spin 
of the ground state for the half filled d shell is S" = 5/2, 
which is what we obtain from the diagonalization of the 
model. The sixfold degeneracy at zero field is broken 
by the combined action of spin-orbit and crystal field. 
Due to the spin-orbit coupling, the total spin S and total 
angular momentum L are no longer good quantum num- 
bers. However, our CI method allows to calculate any of 
these expectations values. We have verified that for our 
calculation for the Mn^+ , (S) « 5/2, while (L) « 0, with 
a deviation smaller than 0.1%. In the same way, Sz is 
almost a good quantum number. 

The location of the first neighbors of Mn is taken from 
referenced. The values of e' and e" are taken so that the 
lowest energy levels of the energy spectra obtained from 
the diagonalization of are in agreement with those 
of the single ion Hamiltonian, as shown in the Figl5l^b). 
At zero field, the lowest energy doublet corresponds to 
(Sz) ~ ±5/2. For the two pairs of excited levels, we get 
(Sz) « ±3/2 and (Sz) ~ ±1/2, in order of increasing 
energy. Fig. [5jb) shows the magnetic field dependence of 
the low energy spectra of the Mn^+ obtained using the 
CI calculation, together with the fitting to a phenomeno- 
logical spin model^Ti^i^iSl 

Hs = DS^z + E{Sl - SD + gfiBB.S. (28) 

The first two terms in Eq. describe the single ion 

magneto-crystalline anisotropy while the last one corre- 
sponds to the Zeeman splitting term under an applied 
magnetic field B. The main magnetization direction z in 
Eq. p8| depends on the substrate and magnetic atom 
nature. In the case of the Mn on a CU2N substrate, the 
z-axis is perpendicular to the surface. This result is also 
reproduced by our model (j26p . 

C. Transport 

Once the many-body eigenstates of Hamiltonian ([25)1 
are obtained, we are in position to study transport 
through the magnetic atom. For that matter, we need 
to specify the coupling of the 5 d orbitals to the tip and 
the substrate. As the dl/dV spectra was recorded with 
the tip located exactly over the Mn atom, we will assume 
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FIG. 5: (Color online) (a) Single particle energy levels of 
the crystal field Hamiltonian with the dominant orbital con- 
tribution (only red levels are coupled by the crystal field). (b) 
Lowest energy spectra of the total Hamiltonian T-Lms for the 
Mn^''' ion on a CU2N surface as studied in Ref. 0- The dots 
symbols corresponds to the solution of the phenomenological 
spin model with D = — 0.39meV and E = O.OBmeV. Here 
Ed = -21.5eV, e' = 11.3 and e" = 1.9. 
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FIG. 6: (Color online) (a) Total (solid line), elastic (dashed 
line) and inelastic (dotted-dashed line) dl/dV as a function 
of applied bias for the Mn^"*" ion {psVs = 1, prVr^sz^^r^ ~ 1, 
T = 0.4K and Ea = -21.5eV). (b) Ratio G'"<=(y)/G=''"'(0) 
for V = 2meV versus on-site energy Ed (lower axis) or (fie — 
Ef) /Ajj. (upper axis), with A/i — fie — fih- Other parameters 
as in Fig. \5\ 



VI. DISCUSSION AND CONCLUSIONS 



that the tip-atom tunneling is dominated by tunneling 
between the tip apex s orbital and the d^z^-r^j^^ ori- 
ented along the adatom-tip axis. For the coupling with 
the substrate, the situation is significantly more com- 
plicated and we couple equally all the d-orbitals to the 
substrate, i.e., Vs,i — Vs- For simplicity, we will omit 
the coupling between the s-orbital of the Cu atom and 
the empty s-orbital of the Mn^"*" ion at odds with exist- 
ing DFT calculation!^ Another important parameter to 
properly account for the transport properties of the sys- 
tem is the Fermi level of the electrodes. Here we have 
assumed that the Fermi level of the Cu substrate coin- 
cides with its bulk Fermi level, Ep = — 7eVi^ 

The resulting dl/dV is plotted in Fig. [SJa), were the 
elastic, inelastic and total differential conductance are 
plotted. Our calculation reproduces both the line-shape 
of the dl/dV curves as well as the the relative contribu- 
tion between the elastic and inelastic parts, G-mei/Gd — 
0.5 4.!^' Within the model this ratio depends on the po- 
sition of the charging energies of the atom, and fj,h, 
with respect to the chemical potential of the electrodes. 
In Fig.injb) we show ratio Ginci/Gd as a function of the 
on-site energy level Ed in the window of energies where 
the system ground state contains 5 electrons. As ob- 
served, the ratio Ginci/Gd varies smoothly between 0.4 
and 0.6. Thus, our model yields a large inelastic signal, 
consistent with the experiments, without fine tuning the 
on-site energy Ed- Notice that in the case of Mn on CU2N 
at T = 0.4K, the thermal broadening of the inelastic step 
is such that the inelastic conductance is non-zero even at 
zero bias. 



We have shown that the inelastic tunneling spec- 
troscopy widely used to study magnetic molecules and 
atoms adsorbed on surfaces, can be understood in terms 
of cotunnelingi^i^ As the electrons go from the tip to the 
surface, the magnetic system must undergo virtual tran- 
sitions to states with an extra electron or an extra hole. 
This picture holds both for elastic tunneling, in which 
case the MS returns to the original state after the vir- 
tual charging process, and inelastic tunnel, for which the 
state before and after the virtual charging are different. 
Thus, the origin of both elastic and inelastic conductance 
is the same, which accounts for the large inelastic signal 
reported experimentally in a variety of systems, includ- 
ing Mn, Fe and Co on Cu2N'^ or Fe on InSbf^^ Further 
support to this claim comes from comparison of the evo- 
lution of the dl / dV as a function of an applied magnetic 
field of a quantum dot with a single resident electron 
in the Coulomb Blockade regimSi^ and a single Cobalt 
atom on CU2N, both undergoing a transition from the 
Kondo regime at low field to inelastic steps at large field. 
Both systems show very similar dl /dV. In addition, our 
microscopic theory provides a natural starting point to 
describe both the appearance of Kondo correlations, and 
their relation to the inelastic spin flips in the context of 
magnetic adatoms and molecules. 

Our approach is based on the derivation of an effec- 
tive cotunneling Hamiltonian acting only in the space 
of neutral configurations of the MS. The calculation of 
the effective Hamiltonian requires the exact diagonaliza- 
tion of the MS in the neutral subspace as well as the 
subspaces with one extra electron and one extra hole. 
From the formal point of view our results are in agree- 
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ment with previous works based on a truncation of the T- 
matrix to second order in the couphng Hamiltoniani^""— 
Our approach permits to obtain an effective cotunnehng 
Hamiltonian that can be compared with effective Kondo- 
hke Hamiltonians proposed in most theoretical analysis 
of lETS experiments^-^-22 

We have also explored the origin of the experimentally 
observed asymmetry with respect to bias inversion in the 
dl/dV curvesi^i^"— It comes from a combination of two 
ingredients. First, we need to consider that the bias volt- 
age results in a shift of the chemical potential in the tip, 
the one in the surface remaining constant. Second, the 
energy level alignment of the MS must be such that one of 
the cotunneling channels, either virtual electron addition 
or virtual hole addition, is dominant. 

In summary, we propose a method to describe sin- 
gle spin inelastic electron tunneling spectroscopy which 
does not rely on effective spin models to describe both 
the magnetic system and the spin-flip assisted tunnel- 
ing. Our approach provides a natural explanation for 
the large inelastic signals observed experimentally, and a 
microscopic mechanism for the spin assisted tunneling. 
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Appendix A: Effective tunneling Hamiltonian 

We now use Eq. Q to derive an effective Hamilto- 
nian which acts on the reservoir fermions and on the 
go subspace of the central island only. By so doing, we 
shall eliminate the and d operators from the effective 
Hamiltonian and, more important, we shall obtain a tun- 
neling Hamiltonian for which the current can be derived 
straightforwardly. The matrix element between any two 
states in the qq manifold can be written as: 

(^|Vtu„|A') = (vl//(0)Kn|Vtun|*/'(0))|n'), (Al) 

where |A^) = \n) (g) |*/(0)), with l^t^/} a multi-electronic 
Slater state describing independent Fermi seas of left and 
right electrodes. Importantly, the unperturbed states are 
product states of the left and right electrodes and the cen- 
tral island. These states can describe both, the ground 
state of the MS with no excitations in the electrodes and 
excited states with an electron-hole pair in the electrodes 
and a excited state n' in the central island. Notice that 
the electron-hole pair can be either in one electrode or 
split in the left and right electrodes. In the second case, 
this excitation contributes to the net current flow. Now 



we need to evaluate matrix elements like 

{^fm{n\V+\M_) = Y.Kd'^fiO)\fa\%nf{-)) 

'x(n|4|m_). (A2) 

Before going further, it is convenient to write down 
the explicit form of the electrodes wavefunctions. If we 
denote the ground state of the electrodes in the Fermi 
sea with no excitations and in its neutral charge state 
as |0), we can write |\l//(0)) = flfa\0), where we are 
creating an electron-hole pair with quantum number a. 
For the states with one electron excess (defect) we will 

have \^r,.fH) - flflUm (|*™/(+)) = fpflfalO)). 
The matrix element of the electrode operator in Eq. (jA2p 
selects one and only one term in electrode part of the 
sums J2m- = Em- Sm,/- term in question is such 
that 



l*™/(-)) = /j|*/(0)). 



(A3) 



This relation is equivalent to write (^'/(0)|/-,.|^m/(— )) — 
(1 - nfij))Sp^, where n^i^)^ = (^/(0)|/t/^|V'/(0))_ is 
the zero temperature occupation of a quasiparticle with 
quantum number 7. We can now write 



M_ 



X „ (^/ (0) I ul \^f, (0)} 



E„ 



X (n|c?i |m_)(r7i_ \di> \n'). 



(A4) 



A similar expression can be obtained for the matrix ele- 
ments involving states 



E. 



Em+ — Eq 

V -V* ■ 

Eq — Cq 



(V'/(0)|/j;/a'|^/'(0)) 



m-i- 



X {n\di\m^) {m^\dl \n'). 



(A5) 



Now, it is straightforward to show that the addition of 
Eqs. (|A4[) - (jA5p leads to the final expression 



Appendix B: Tunneling transition rates 



As stated in the Sec. Ill Bl the cotunneling transition 
rates can be calculated applying the Fermi Golden Rule 
to the effective tunneling Hamiltonian T-Lcotun ■ Introduc- 
ing the density of states p,ja and using Eqs. (|A4IIA5j) . the 
transition rate from a state n of the central island to an 
state n' , with the transport electron going from electrode 
77 to 77' and its spin from a to a', are given by 
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2tt f 

,(6,e + A)|n') - {n'\d';-},(e + A,e)\n) ^ 



(Bl) 



with Ann' = En — En' and /(e) the Fermi-Dirac distribu- 
tion. The matrix elements of the O'^^'^ operators in Eq. 
(IBll) are defined as 



("ic>i:.U'(^'^')i-')= 



E 

ii' ,m.j^ 



En 



Eo 



nil r ■ -I l\ 



where we have used a simpUfied notation Vrf,i(t) = 
^fc(e)r).i -These transitions rates are in perfect agreement 
with the rates obtained by a second order truncation of 
the T-matrix.i^^ 



and 



E 



yu^)yr,'A^) 



Em_ — En 
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